/* All the quantities are defined in SI-units unless specifically stated. */

#ifndef SCENE_H

# define   SCENE_H
# define   H_FILE  "scene.h"

# include  <stdio.h>

/* Overall, parallelization can be controlled using this parameter. */
# ifdef _OPENMP
// #    define   _OPENMP_
# endif

/* Paralallization of a function can be controlled using these parameters. */
# ifdef _OPENMP_
#    define   _OMP_VEL_ADVANCE_
#    define   _OMP_VEL_BC_
#    define   _OMP_SOLVE_POISSON_
#    define   _OMP_PRESS_BC_
#    define   _OMP_TEMP_ADVANCE_
#    define   _OMP_TEMP_BC_
# endif

/* This flag sets what version of the solve_poisson should be called. If it is defines, OpenMP version is called. */
// # define CALL_SOLVE_POISSON_OMP

# define   SOLVE_ENERGY

/* If this parameter is defined, various configuration parameters need to be confirmed by user in the
   begining of the program. */
// # define   CHECK_CONFIG

/* Maximum filename length. */
#  define   MAX_FNAME_LEN    50

/* It should be zero if the program statrs from begining. It can be 0, 1, 2,....9. */
#  define   RUN_CNT  0

/* This specifies if the initial condition should be read from a recovery-file or not. */
// #  define   READ_RECOV_FILE

/* Recovery file: data may be read from this file as initial condition for velocity, temperature and pressure fields.*/
#  define   RECOV_FILE          "recov"

/* If poisson equation takes more than this number of iteration for convergence, a recovery file is written. */
#  define   RECOVF_POISS_CNT    100000

#  ifdef READ_RECOV_FILE
 /* 200 the maximum length of a line in the 'RECOV_FILE' */
#   define MAX_RECOVF_LINELEN  200
#  else
 /* The initial time-offset when calculation starts (in second). */
#   define t_INIT    0.0

 /* The constant time-step for time marching (in Seconds). */
#   define DELTAt    2.0E-5
#  endif

/* The maximum time until which the calculation has to continue (in seconds). */
# define tMAX        7200.0

/* After the time tFRAC, the flow rate reaches to wFRAC-th of the maximum possible flow-rate.
   For example, if Q(t = 0) = 0; and, Q(t = infinity) = Qmax;  then, Q(t = tFRAC) = wFRAC*Qmax.
   This also implies that: W_center(t = tFRAC) = wFRAC*W_center(t = t_infinity). */
# define tFRAC  1.0    /* In Seconds. */
# define wFRAC  0.99

//# define   r_UNIF_GRID
# define   th_UNIF_GRID
//# define   z_UNIF_GRID

/* Files in which grid distances are stored. */
# define   R_FILE   "r-grids"
# define   TH_FILE  "th-grids"
# define   Z_FILE   "z-grids"

/* If this variable is defined, the grids are read from the file R_FILE, TH_FILE and Z_FILE. */
// # define   READ_GRID_FILES

/* Various files used to record statistics and informations during program run.
   LOG_FILE     :  To store essential information like number of grids, grid-size etc.
   ittr_RES_FILE : To store different iteration-parameters at every time-step.

   If any of these files are defined to (char *)0, the output is set to "stdout". If it is not defined,
   nothing will be printed to the file. */
# define  LOG_FILE         "log"      /* "log",   (char *)0 */
# define  ittr_RES_FILE    (char *)0     /* "ittr",  (char *)0 */
/* ittr_RES_FILE will be renewed after IRFILE_FREQ itterations. */
# define  IRFILE_FREQ       200000
/* After ittr_ENTRY_FREQ, there will be a new entry in the ittr_RES_FILE. */
# define  ittr_ENTRY_FREQ   100
/* Until first ittr_PRINTALL_CNT, the parameter ittr_ENTRY_FREQ will be ignored and entries will be made at every time-step. */
# define  ittr_PRINTALL_CNT 200
/* The frequency with which the pressure-residuals will be printed in iit_RES_FILE while solving the PPE. */
# define  PPE_RESUAL_FREQ   100

# define  RES_FILE   "res" /* The result files. */
# define  RESF_FREQ   5000

/* Number of dependent variables. Do not edit it. */
# define   DEPVAR_N  4 /* u, v, w, P */
# ifdef    SOLVE_ENERGY
# undef    DEPVAR_N
# define   DEPVAR_N  5 /* u, v, w, P, T */
# endif

/* Runge-Kutta related parameters. Do not edit it arbitariry. */
# define   nRK   3    /* No. of stages and order of the RK method. */
# define   c2    (1.0/6)  /* The free parameter used to determine the co-efficients. */
# define   a31   (2.0/3 - 2.0/(9*c2))
# define   a32   (2.0/(9*c2))
# define   a41   0.25
# define   a42   0
# define   a43   0.75
# define   b1    a41
# define   b2    a42
# define   b3    a43
# define   d1    (-2+1.0/c2)
# define   d2    (-1.0/c2)
# define   d3    3.0


/* PPE relaxation factor */
# define   PPE_RFACT   1.65

/* If pressure-residual during solution of the Poisson's equation in 'solve_poisson()' remains less than 'PRESS_DELTA'
   for 'POISS_CNT_MIN' times continuously, we say that the Poisson Equation is solved. This also implies that the
   "p_cnt" returned by 'solve_poisson()' will never be less than 'POISS_CNT_MIN'.

   "POISS_CNT_MIN >= max(nr, nth, nz)" should be used.

       If the above condition doesn't meet within 'POISS_CNT_MAX' iterations, we say that Poisson Equation
   could not be solved as expected. The details can be found in 'T_RES_FILE'. */
# define   PRESS_DELTA      1.0E-4
# define   POISS_CNT_MIN    10
# define   POISS_CNT_MAX    200000

/* "solve_poisson_omp()" ignores my this flag because this function always uses JACOBI itteration method.\
    However, "solve_poisson()" can be controlled using it. */
// # define  PPE_ITTR_JACOBI


/* If the initial-residual during solution of the Poisson Equation is less than "PRESS_DELTA" (i.e. "solve_poisson()"
   returns "POISS_STOP_CNT") for "SST_CNT_MAX" times continuously, the CFL-number is increased.
        If Poisson Equation doesn't converge within "POISS_CNT_MAX" for "DIV_CNT_MAX" times continuously, we say that
   the solution is diverging and try to get it converged by decreasing the CLF-number.

   The aforementioned increament or decreament in CFL-number is done only if "CLF_MIN <= cfl <= CLF_MAX". If this
   condition doesn't hold true, the solution is assummed to reached steady state or diverged respectively and the
   program is stopped.

   To disable usage of these, set these to a negative number.
*/
# define   SST_CNT_MAX      -1
# define   DIV_CNT_MAX      1

/*
# ifndef  DELTAt
#   define   CFL_INIT      0.02
#   define   CFL_MIN       0.01
#   define   CFL_MAX       0.1
#   define   CFL_DELTA     0.001
# endif*/


# ifdef SOLVE_ENERGY
 /* This factor controls the dissipation added to convective terms of
    the energy equation.
    If EEQ_DISSIP_FACT = 1, it corresponds to full upwinding scheme.
    If it is zero, it corresponds to full central difference scheme.
    If it is undefined, its value is considered as '1' and the convective
    term is discretized in simplified way. */
#   define  EEQ_DISSIP_FACT  0.5


#   define  ABS_ZERO     273.0  /* in Kelvin. */
#   define  T_INIT       (21.5 + ABS_ZERO)  /* Initial temperature in Kelvin. */
#   define  T_INLET      (21.5 + ABS_ZERO)  /* Inlet temperature in Kelvin. */


 /* If VOLUME_HSRC = 'y' is defined, the heat-source-term will be added to the energy-equation, otherwise, not. */
#   define  VHSRC_FLAG       'n'
#   define  EXPECTED_T_RISE  5.0  /* In Kelvin */
#   define  EXPECTED_TIME    600.0  /* In Sec */


/* All the solid surfaces are considered to be insulating excluding the central region of the top-surface where
   a constant heat-flux heater is attached. All the heat-fluxes are in "W/m^2". The heater is assummed to be on
   after "t_HEATER_ON" seconds from begining of the flow. */
#   define  t_HEATER_ON      0.0
#   define  TOP_HEATER_HF    20.0E+3
#   define  TOP_ANNULUS_HF   0.0
#   define  SIDE_WALL_HF     0.0
#   define  SOLID_BASE_HF    0.0

 /* These three parameters defines if Kinematic Viscosity, Density and Thermal Conductivity of water are temperature-dependent or not. */
#   define  KVISC_TDEP
#   define  DENS_TDEP
#   define  THMCOND_TDEP

/* If RETURN_CONST_KVISC and RETURN_CONST_DENS are defined, the functions kvisc() and dens() will return
   values at T = T_INIT. *//*
#   ifdef  KVISC_TDEP
#      define RETURN_CONST_KVISC
#   endif

#   ifdef  DENS_TDEP
#      define RETURN_CONST_DENS
#   endif

#   ifdef  THMCOND_TDEP
#     define RETURN_CONST_THMCOND
#   endif */
# endif

# ifndef  KVISC_TDEP
#   define   KIN_VISC      0.999028E-6  /* In "m^2/sec". */
# endif
# ifndef  DENS_TDEP
#   define   DENSITY       997.97       /* In "Kg/m^3". */
# endif
# ifndef  THMCOND_TDEP
#   define   COND_WATER    0.6          /* Water conductivity in 'W/m.K'. */
# endif
# define   Cp_WATER      4.172505E+3  /* Specific heat of water in 'J/Kg.K'. */
# define   W_MAX         0.12     /*0.121709*/     /* For parabolic profile at the inlet, Re = 600. */
# define   GRAVITY       9.81         /* Gravitational constant. */
# define   P_ATM         1.01E+1      /* In "Pa". */


/* Correct w_outlet for initial N_w_CORRECT number of steps so that "dmdt_in = dmdt_out". */
// #  define N_w_CORRECT 3

// # define   CORRECT_T_OUTLET

/*  If tn_BACKUP = 3, it stores physical quantities at time-instants of "t - dt" and "t - 2*dt in
    recovery file. Also, quantities at t, t-dt, and t-2*dt are used in application of Orlansky BC. */
# define  tn_BACKUP      3

/* Physical-domain velocities are stored specifically at a different memory location for zCells_BACKUP
   number of starting from the cell adjascent to the outlet boundary continuing along z-axis inside
   the physical domain. If zCells_BACKUP = 3, the quantities are saved for cell numbers NzG, NzG+1
   and NzG+2 at the outlet boundary. */
# define  zCells_BACKUP 3

/* This defined the number of sections the whole domain is divided into for the computation. */
# define Nsec  3

struct arr_struct
   {
    double ****q;
    unsigned short nq, n1, n2, n3;
   };

struct orlansky
   {
    double ***u[tn_BACKUP], ***v[tn_BACKUP], ***w[tn_BACKUP], **tnext_ug, **tnext_vg, **tnext_wg;
#   ifdef SOLVE_ENERGY
    double ***T[tn_BACKUP], **tnext_Tg;
#   endif
   };

# define  N0rZONE  2
# define  N1rZONE  6
# define  N2rZONE  1
# define  NrZONE_MAX  N1rZONE

# define  N0zZONE  3
# define  N1zZONE  3
# define  N2zZONE  4
# define  NzZONE_MAX  N2zZONE

struct grid_config
 {
# ifndef r_UNIF_GRID
  double r_cd[NrZONE_MAX], dr[NrZONE_MAX+1], r[NrZONE_MAX+1];
  short nr[NrZONE_MAX], nzone_r;
# endif

# ifndef z_UNIF_GRID
  double h_cd[NzZONE_MAX], dh[NzZONE_MAX+1], h[NzZONE_MAX+1];
  short nz[NzZONE_MAX], nzone_z;
# endif
 };


struct grid
 {
  unsigned short nr_nozzle, nr_gap, nr[Nsec], nth, nz[Nsec], nTRD;
  double *r[Nsec], *th, *z[Nsec], ***p[Nsec], ***u[Nsec], ***v[Nsec], ***w[Nsec], ***um[Nsec], ***vm[Nsec], ***wm[Nsec];

# ifdef r_UNIF_GRID
  double dr;
# endif

# ifdef th_UNIF_GRID
  double dth;
# endif

# ifdef z_UNIF_GRID
  double dz;
# endif

# ifdef SOLVE_ENERGY
  unsigned short nr_heater;
  double r_heater, ***T[Nsec];
  char vhsrc_flag;
# endif

  double dt, t, density0, kvisc0, Re0, thmcond0;
  FILE *lptr, *trptr, *tptr;

# ifndef DELTAt
  double cfl;
# endif

  struct orlansky outBC;
 };


/* Number of ghost-cells. */
# define     NrG        2
# define     NthG       2
# define     NzG        2

/* Degree of the polinomial used to approximate the ghost-cell w-velocity
   at the outlet. Used in "vel_bc()". */
// # define     VEL_INTRP_N    2

/* Degree of polynomial used to approximate the ghost-cell pressure
   at the outlet. Used in "press_bc()". */
# define     PRESS_INTRP_N  2

# ifdef SOLVE_ENERGY
 /* Degree of polynomial used to approximate the ghost-cell temperature at inlet in "temp_bc()". */
#   define   TEMP_INTRP_N   2

/* DEFINING r-COORDINATE RELATED QUANTITIES. */
#   define   d_HEATER   100.0E-3 /* "0 <= d_HEATER <= d1_CYL"; Modified. */
#   define   r_HEATER   (0.5*d_HEATER) /* Modified. */
# endif
# define     d_NOZZLE   9.85E-3
# define     d1_OUTLET  233.5E-3 /* Modified. */
# define     d1_CYL     273.5E-3 /* Modified. */
# define     d0_CYL     40.0E-3
# define     r_NOZZLE   (0.5*d_NOZZLE)
# define     r1_OUTLET  (0.5*d1_OUTLET) /* Modified. */
# define     r0_CYL     (0.5*d0_CYL)
# define     r1_CYL     (0.5*d1_CYL) /* Modified. */

# ifdef  r_UNIF_GRID
#   define   dR_CYL     0.985E-3 /* Modified. */
# else
#   define   dr_CENTER  0.3E-3
#   define   dr_NOZZLE  0.4E-3 /* Modified. */
#   define   dr0_CYL    0.6E-3 /* Modified. */
#   define   r1_MID1    22.0E-3  /* Modified. */
#   define   r1_MID2    90.0E-3  /* Modified. */
#   define   dr1_OUTLET 0.8E-3
#   define   dr1_CYL    1.3E-3 /* Modified. */
# endif


/* DEFINING theta-COORDINATE RELATED QUANTITIES. */
# ifdef th_UNIF_GRID
#    define  dTHETA     5.0  /* modified */
#    define  nTH_COMP   5 /* It should be more than 2*NthG. */
# endif


/* DEFINING z-COORDINATE RELATED QUANTITIES. */
#ifdef z_UNIF_GRID
#    define  dh            0.975E-3 /* Modified */
#else
/* For the top-section 1. */
#   define   h1_CYL        18.6E-3
#   define   dh1_TOP       0.1E-3
#   define   dh1_MID2      0.5E-3 /* Modified. */
#   define   h1_MID2       12.0E-3
#   define   h1_MID1       8.0E-3
#   define   dh1_MID1      0.8E-3 /* Modified. */
#   define   dh1_BOT       0.5E-3 /* Modified.  */

/* For the bottom-section, i.e. SECTION-0. */
#   define   h0_CYL        40.0E-3
#   define   h0_MID2       30.0E-3 /* Modified */
#   define   h0_MID1       10.0E-3
#   define   dh0_MID1      1.0E-3 /* Modified */
#   define   dh0_BOT       0.5E-3 /* Modified */
  /* Height of the water level in the tank from the nozzle-exit. */
#   define   h_WATER       (h0_CYL+h1_CYL)

/* For the extended section-2 at the outlet to take care of the back-flow. */
#   define   h2_CYL        20.0E-3
#   define   h2_MID2       15.0E-3 /* Modified. */
#   define   h2_MID1       8.0E-3 /* Modified. */
#   define   dh2_MID1      1.1E-3 /* Modified. */
#   define   dh2_UNIF_BOT  0.6E-3 /* Modified. */
#   define   nz2_UNIF_BOT  1
#endif

#endif
